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A Segmentation  Method  under  Geometric 
Constraints  after  Pre-processing 

D.  Apprato,  J.  B.  Betbeder,  C.  Gout,  and  S.  Vieira- Teste 


Abstract.  For  a geophysical  image  with  homogeneous  grey  levels,  we 
propose  a method  of  segmentation  that  could  be  subdivided  into  two  parts: 
the  first  one  concerns  a pre-processing  of  the  image  which  provides  an 
enhancement  of  some  features  present  on  the  image.  The  originality  of  the 
method  consists  in  using  a scale  transformation  applied  to  the  pixel  values 
of  the  image.  The  second  part  presents  a segmentation  method  using 
deformable  surfaces.  The  originality  of  this  segmentation  method  is  that 
it  considers  the  active  contour  model  as  a set  of  articulated  curves,  which 
corresponds  to  the  interfaces  between  different  layers  and  faults.  Moreover, 
the  a priori  knowledge  of  well  data  allows  us  to  make  some  geometric 
constraints  on  the  model.  The  solution  is  obtained  by  minimization  of  a 
nonlinear  functional  under  constraints  in  a suitable  convex  set.  Solving  the 
minimization  problem  consists  in  particular  in  a fc-order  Taylor  formula 
applied  to  linearize  the  nonlinear  term. 


§1.  Segmentation  Pre-processing 

Image  segmentation  is  one  of  the  most  important  steps  leading  to  the  analysis 
of  processed  image  data.  Its  main  goal  is  to  divide  an  image  into  parts  that 
have  strong  correlation  with  objects  or  areas  of  the  real  world  contained  in  the 
image.  The  image  is  divided  into  separate  regions  that  are  homogeneous  with 
respect  to  a chosen  property  such  as  brightness,  color,  reflectivity,  context,  etc. 
However,  in  certain  cases,  the  grey  levels  of  an  image  could  be  homogeneous 
and  make  the  segmentation  more  difficult  to  realize.  This  is  particularly  true 
in  the  case  of  geophysical  and  medical  images  (cf.  [14,15]).  In  the  first  part  of 
this  work,  we  propose  a method  to  solve  this  problem  using  families  of  scale 
transformations.  The  use  of  scale  transformations  is  common  in  imaging.  The 
aim  of  this  pre-processing  is  an  improvement  of  the  image  function  data  that 
suppresses  unwilling  distortions,  or  enhances  some  image  features  important 
for  further  processing.  It  provides  improvement  of  the  contrasts,  and  it  rep- 
resents a tool  to  pre-process  images  used  in  most  computer  algorithms  today. 
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According  to  Sonka,  Haclav  and  Boyle  [15],  the  pre-processing  of  images  may 
be  classified  into  four  categories  (pixel  brightness  transformations,  geometric 
transformations,  pre-processing  methods  that  use  a local  neighborhood  of  the 
processed  pixel,  and  image  restoration  that  requires  knowledge  about  the  en- 
tire image)  according  to  the  size  of  the  pixel  neighborhood  that  is  used  for  the 
calculation  of  a new  pixel  brightness.  The  transformation  of  the  brightness 
and  of  the  contrast  of  an  image  allows  us  to  focus  on  phenomena  that  are  hard 
to  see  in  the  plain  image. 

For  a given  image,  we  are  going  to  consider  the  pixel  values  as  a topo- 
graphic map:  the  brightness  value  of  each  pixel  is  the  height  of  the  (hyper-) 
surface  at  that  point.  For  a data  set  of  pixels  (xj,y,-,  Zi,  A (xj,t/;,  Zi)){,  we 
apply  the  following  functions: 

. C d-  A{xi,yi,zi)i  C [0,255]  — ♦ [0,255], 

• r'(woM))eff”(fi,R), 

• ipd  (Td  ((fid  o (qdoA)))  converges  to  C o A when  d converges  to  0, 

where  A is  an  attribute  function  introduced  in  Section  2.1,  qd  (resp.  ipd  and 
il>d)  are  scale  transformations  converging  to  c (resp.  ip  and  ip),  and  Td  is  a 
Dmspline  operator  (see  Arcangeli  [2]). 

The  scale  transformation  qd  converges  to  a usual  brightness  transforma- 
tions ? (see  Apprato  and  Gout  [1]):  for  instance,  C could  be  a scale  transfor- 
mation which  enhances  the  image  contrast  between  brightness  values  pi  and 
P2- 

Let  us  consider  the  subdivision  {■u1,U2,  ...,«*,  •••,  «P(d)  }i=i  of  the 

interval  [0,255]  satisfying  C (A  (x*,?/,-,  2j))  = Ui,  p(d)  being  the  number  of 
different  pixel  values  of  the  image  (<  255  for  a grey  scale  image).  The  function 
Cd  is  defined,  for  any  x € [A(xi,yi,Zi)  = Wi,A(xi+1,yl+uzi+ 1)  = wi+i\,  and 
for  an  integer  1 < i < p{d)  — 1,  by 

Cd(x)  = Uiq$m  [(x  - Wi)  / (wi+1  - Wj)]  + ui+1qlm  [(x  - W;)  / (wi+1  - w*)] 

+ Qi  (wi ) (wi+ 1 - Wi)  q°lm  [{x  - Wi)  / (wi+i  - iOj)] 

+ a i {wi+ 1)  (wi+i  - Wi)  q\m  [(x  - wt)  / (wi+1  - w;)] , 

where  the  q\m,  for  l = 0, 1,  and  j = 1,  ...,m,  are  the  Hermite  finite  element 
basis  functions,  and  where  ai  (w*)  = (u;+ 1 — w,)/  (wj+i  — w,)  and  aq  ( wp (d))  = 
(wp(d)  ~ wp(d)- 1)  / (Mp(d)  _ up(d)~ l)-  Then,  Gout  [9]  showed  that  for  any  d € 
D,  for  an  integer  i,  1 < i < p(d)  — 1,  Cd  (w;)  = and  Q € Cm([0,  255]). 

Likewise,  in  order  to  recover  a finer  image,  it  is  useful  to  apply  the  “large 
variations”  algorithm  introduced  in  [9].  In  fact,  after  having  applied  the  func- 
tion Cd  to  improve  the  contrast  of  the  image  (and  thus  increasing  the  variations 
of  the  corresponding  pixel  values),  it  is  very  useful  to  use  a method  that  takes 
into  account  these  rapidly  varying  data.  Let  us  note  that  even  without  using 
the  scale  transformations  Cd,  an  image  often  has  large  variations  (this  occurs 
for  example  when  a dark  zone  is  near  a brighter  one).  That  is  why  we  propose 
to  use  the  “Large  variations”  algorithm.  This  algorithm  uses  two-scale  trans- 
formations, namely  ipd  for  the  pre-processing,  and  ipd  for  the  post-processing. 
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The  first  one,  <pd,  is  used  to  suppress  the  oscillations  of  the  data.  The  pre- 
processing function  tpj.  is  such  that  the  data  do  not  present  large  variations, 
and  therefore  a usual  spline  operator  Td  (e.g.  [2])  can  subsequently  be  applied 
without  generating  significant  oscillations.  The  second  scale  transformation 
ipd  is  then  applied  to  the  approximated  values  to  map  them  back  and  obtain 
the  initial  approximated  pixel  values.  It  is  important  to  underline  that  the 
proposed  scale  transformations  do  not  create  parasitic  oscillations.  Moreover, 
this  method  is  applied  without  any  particular  knowledge  of  the  location  of  the 
large  variations  in  the  dataset. 

So,  for  pre-processing,  we  propose  two  algorithms:  in  the  first  one,  we  just 
apply  a scale  transformations  Q as  a brightness  transformation  for  contrast 
enhancement,  in  the  second  one,  we  also  apply  the  “large  variations”  algorithm 
in  order  to  obtain  a finer  represention  of  the  image  which  represents  the  main 
advantage  of  this  approach. 

The  reader  is  referred  to  [1,8,10]  for  a complete  study  of  this  method, 
including  its  convergence  and  numerical  results.  Let  us  note  that  this  method 
is  also  efficient  for  noise  removal  as  shown  in  [1]. 

§2.  Segmentation  Method 

We  use  deformable  models  (external  forces,  evolution  term,  see  Kass,  Witkin 
and  Terzopoulos  [16,17])  and  classical  approximation  techniques  such  as  spline 
theory  (see  de  Boor  [3],  Laurent  [12],  Schumaker  [13])  and  the  finite  element 
method  [4]. 

We  propose  an  analytic  approach  which  uses  deformable  models  instead 
of  a geometrical  one  as  done  for  instance  in  Sethian  [14].  We  recall  that  the 
principle  of  the  deformable  model  method  lies  in  attracting  the  representation 
towards  the  structure  using  forces: 

• Internal  forces  describing  properties  of  elasticity  and  rigidity  of  the  rep- 
resentation, connected  to  its  derivatives  (e.g.,  the  energy  of  thin  plates); 

• External  forces  coming  from  potentials  which  characterize  the  elements 

of  the  structure  with  respect  to  the  attributes  data. 

Geometrical  constraints  are  associated  with  well  interpolation  conditions 
(case  of  geophysical  images  with  well  data).  Deformable  models  provide  a way 
of  interactively  acting  on  the  representation  by  adding  a dynamic  term  in  the 
minimization  problem  (see  for  instance  Cohen  and  Cohen  [5],  Cohen,  Cohen 
and  Ayache  [6],  and  Cohen,  Bardinet  and  Ayache  [7]),  that  permits  upgrading 
the  models  to  the  solution  of  the  minimization  problem  introduced. 

In  this  section,  we  first  give  the  geophysical  data  and  then  the  minimiza- 
tion problem  is  studied.  The  nonlinear  problem  and  its  discretization  are 
given  in  the  subsequent  sections. 

2.1.  The  data  on  the  structure 

Two  types  of  data  are  available:  attribute  data  and  well  data.  For  each  at- 
tribute A,  the  attribute  data  are  (xi,yi,Zi,A(xi,yi,Zi)),  where  {xi,yi,z{)  are 
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the  coordinates  of  the  barycentre  of  a voxel,  and  A z{)  is  the  attribute 

value  A in  this  voxel.  The  well  data  are  depth  data: 

(Xj,yj,zj)j=1  N = aj  where  N is  the  number  of  interpolation  points.  This 
model  allows  a conceptual  representation  of  the  structure  by  identification 
of  its  various  elements,  and  permits  topological  connections  between  those 
elements.  This  model  induces  the  parameterization  of  the  structure.  Each 
element  of  the  structure  (layer,  fault,  etc.)  is  identified  with  a connection  of 
four  points  with  a label  E.  Furthermore,  each  quadruplet  is  connected  by  two 
points  (which  can  be  thought  as  a common  side  of  a “quadrilateral”  repre- 
sented by  the  four  points)  with  another  quadruplet.  Practically,  the  a priori 
model  can  be  constructed  by  introducing  a 3D  block  and  a regular  grid  of  this 
block.  The  aim  is  to  find  a space  of  admissible  representations  consistent  with 
the  a priori  model  and  the  criteria  connected  with  the  data.  Therefore,  it  is 
necessary  to  choose  a space  of  functions  characterized  by  a domain  of  defini- 
tion connected  with  the  a priori  model  and  regularity  conditions  connected 
to  the  data.  The  idea  is  to  transform  the  a priori  model  into  a normalized 
model  called  the  model  of  reference  (denoted  by  M').  For  example,  we  can 
choose  M'  c ft  = [0, 1]  x [0, 1]  x [0, 1].  The  model  M'  is  then  the  image  by 
transformations  of  the  set  of  vertical  and  horizontal  closed  sides  of  the  a priori 
model  as  done  in  [18].  Let  7 be  the  union  of  the  common  edges  of  any  two 
sides  of  M\  we  define  by  M the  interior  of  M'  \ 7.  All  the  functional  spaces 
needed  in  this  work  are  given  in  Vieira-Teste  [18] . 

2.2.  Minimization  criterion 

2.2.1.  Internal  forces:  The  criterion  associated  with  the  internal  forces 
is  a classical  one.  Modelling  this  criterion  bring  us  to  the  following  energy 
functional:  for  any  v € V = H2  (M,  1R3)  fl  C°  (M',IR3), 

Ei  {v)  = M + [v\l  M , 


where 


with  £i  (E)  >0,  Vi  = 1,2,  VE  € M.  The  term  [njj  M corresponds  to  an 
approximation  of  the  elastic  deformation  of  the  model  while  the  term  [i’]2  M 
corresponds  to  an  approximation  of  the  rigid  deformation  of  the  model  (cf. 
Cohen,  Cohen  and  Ayache  [6]). 
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2.2.2.  External  forces:  External  forces  are  issued  from  potentials  connected 
with  attributes.  We  introduce  the  following  energy,  for  any 

£2  (v)  = Y f Px(v/s(s,r))dsdr, 


where  Ps  is  the  potential  associated  with  the  element  parameterized  by  E. 

The  modelling  we  propose  consists  in  minimizing  the  previous  energies  E\ 
and  E2  as  we  will  see  in  subsection  2.2.4.  In  the  case  of  the  velocity  attribute, 
we  use  the  following  potential  to  define  the  layers: 


P(x,y,z)  = -a 


Vl(x,y,z) 


a > 0 


where  A is  the  attribute  “velocity  of  propagation  of  the  seismic  wave” . 


2.2.3.  Interpolation  data:  If  we  suppose  some  parameterization  ( Sj,rj ) € 
M of  each  interpolating  point  aj  = ( Xj,yj , zf)  is  known,  then  we  require  that 

satisfies  v ( Sj,rj ) = aj  for  any  j = 1, ...,  N . 

2.2.4.  Minimization  criterion:  Using  the  notation  and  definitions  intro- 
duced above,  we  consider  the  functional  E defined  on  V by 

E(v ) = Mi ,m  + M 2,M  + Y / ( v/z(s,r ))  dsdr 

ZCMJs 

for  any  v € V.  We  consider  the  set  K associated  with  the  interpolation 
constraints,  and  defined  by 

K = {veV,  Vj  = 1, ...,  N , v(sj,rj)  = aj}. 

This  set  is  convex  and  closed  in  V.  We  also  introduce  the  following  linear 
mapping  (continuous  on  V with  the  norm  j|-||2  M) 

po  : v € V p0v  = (v(8j,rj))j=1  ^ € (IR3)^  . 

We  consider  the  following  minimization  problem:  find  a £ K satisfying 

MveK,  E{a)<E{v). 

We  note  that  this  problem  is  nonlinear  on  the  convex  set  K with  respect 
to  <7.  There  are  two  techniques  to  treat  this  problem.  The  first  one  consists  in 
linearizing  the  nonlinear  term  (linked  to  the  potentials)  in  the  functional  E. 
The  second  one  consists  in  using  the  deformable  models  technique  as  done  in 
the  following  subsection:  we  suppose  that  the  solution  is  a function  of  time, 
which  leads  to  a new  evolution  problem  that  will  be  discretized  both  in  time 
and  space. 
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2.3.  The  nonlinear  problem 

In  this  subsection,  we  give  the  nonlinear  minimization  problem  and  its  dis- 
cretization. Let  us  recall  that  the  deformable  models  technique  consists  in 
assuming  that  a depends  on  time,  and  so  consists  in  adding  a dynamic  term 
to  the  functional  E (<r) 

\§~t  Jm£  (M) (*’ s’ r)  dsdr> 

where  (e  (M)) ^ — e (£)  > 0.  This  term  allows  the  control  at  each  time  of  the 
deformation  of  the  surfaces. 

2.3.1.  Evolution  problem:  Let  T > 0.  We  note 

bL(0,T,  V)  = |w  G L2  (]0,T] , V) , G L2  (]0,T],V')}. 

We  then  consider  the  following  evolution  problem  defined  on  [0,T],  For  any 
t E ]0,T]  and  any  u>  E W(0,T,V),  find  o E W(0,T,  V),o  (t)  € K,  satisfying 

(Pt): 

E^+\§t  / £ ^ <j2 s’ r ) dsdr  - E M + \^t  J f£(M)w2  (*> r ) dsdr< 

with 

o (0)  = cr0  E L2  (M,H3). 

We  are  currently  studying  existence  and  uniqueness  of  (Pt)  using  a Lipschitz 
approximation  of  the  sign  function. 

Likewise,  for  any  t E ]0,  T] , we  consider  the  term 

Lc r(t)  W = ~ (u/s(s>r))  dsdr- 

EC M ,'S 

The  variational  formulation  of  the  problem  (Pt)  with  Kuhn  and  Tucker’s 
relation  is,  taking  as  test  function  v on  the  stationary  space  V (necessary 
condition  without  uniqueness),  for  any  t E ]0,T]  and  any  v E V,  find  (o',  A)  E 
W{0,T,V)  X C°  ([0,T],(]R3)N)  , cr  (t)  E K,  satisfying  (P): 

e (M)  9a^’  —v{s,  r)dsdr  + a(a(t),v)  + (A  (t) , p0v) N<3  = La{t)  (v) 
under  conditions 

ff  (0)  = o0  G L2  (1R3) 

A (0)  = A0  G (H3)", 


and 
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where 

-a(u,t>)  = [»£„  + [«] \M. 

2.3.2.  Discretization  in  time:  In  this  subsection,  we  discretize  (P)  both 
in  time  and  space.  The  originality  of  this  discretization  consists  in  using  a 
k-order  Taylor  development  which  allows  us  to  take  into  account  many  more 
voxels  and  so  to  improve  the  convergence  of  the  method  (see  Vieira- Teste  [18] 
for  more  details).  We  cut  the  interval  ]0,T]  into  sub  intervals  with  length  At. 
Consider 

tm=mAt,  m = \,...,Dx- 

We  use  the  following  approximation  of  the  time  derivative: 

da  u , a{tm)  - <r(tm-i) 
dt  [m)  ~ At 

Assuming  that  am  = a(tm)  and  A™  = A (tm) , we  approximate  the 
variational  problem  as  follows:  For  any  m = 1, Dt  and  any  v £ V,  find 
(<x™,  Am)  e V x (IR3)^  ,cr™  € K,  satisfying  (Pm): 

J e(M)amvdsdr  + At  Ja(<rm,  v)  + (Am,p0n)N  3j 
= f e (M)  am~1vdsdr  + AtL^m  (v) 

Jm 

with  cr°  = <70  € L 2 (M,1R3)  and  A0  = A0  e (H3)N . 

The  previous  problem  is  implicit  and  nonlinear  with  respect  to  the  so- 
lution crm.  We  propose  to  replace  La™  (v)  = La^v  (tTO)  by  a Taylor  series 
expansion  of  order  k > 0 about  the  time  tm.  We  suppose  that  cr  is  in 
Ck  ([0,  T] , L 2 (M,  IR3)) . We  have 

Lam  (^)  — IJa,v  {t-rn) 

and  La<v  ( tm ) ~ La^v  (tm_i)  + AtDLa,v  (tm_i)  + ^^-D2Lc,tV  (tm_i)  + • • ■ + 

Dk La  v (tm_i).  We  note  that  the  problem  (Pm)  is  linear  and  explicit 
with  respect  to  om.  The  following  result  is  based  on  the  Lax-Milgram  Lemma. 

Theorem.  The  problem  (Pm)  has  a unique  solution  (cr™,  A™) . 

2.3.3.  Discretization  in  space:  Let  H be  a nonempty  bounded  subset  in 
1R+  for  which  0 is  an  accumulation  point.  For  any  h £ H,  we  solve  the  min- 
imization problem  (Pm)  in  the  finite  element  space  (14) 3 C V.  The  generic 
finite  element  are  the  Hermite  finite  element  of  class  C 1 for  snakes  and  the 
Bogner-Fox-Schmit  finite  element  rectangle  of  class  C 1 (see  [4])  for  deformable 
surfaces.  To  have  (V/,)3  C V,  it  is  necessary  to  have  a C°  connection  on  7.  To 
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do  that,  it  is  sufficient  to  divide  some  degrees  of  liberty  connected  to  deriva- 
tives as  done  in  [18]. 

We  denote  by  (a™ , ) the  coordinates  of  <rm  in  the  basis  of  Vh  and 
by  (A™, AJy)  the  coordinates  of  Am  ( Mh  — dim(Vh)).  If  <rm  is  a solution 
of  the  discretized  problem  ( Pm ) in  (Vh)3,  we  can  write  am  in  the  basis  of 
(Vh)3:  Vm=l,..,DT,  Vg  = 1,2,3, 


Mh 

(*"T=x>r)V- 

j=i 

with  (a™)'1  € H and  where  the  (<Pj)j=\  Mh  are  t^le  t>asis  functions  of  Vh- 
In  the  following,  we  miss  out  q and  h.  Taking  v = tpi  in  (P),  we  have 
to  solve  (for  q = 1,2,3,  in  the  linear  problem  (P))  a system  of  (Mh  + N) 
equations  with  (Mh  + N)  unknowns.  We  easily  show  that  this  system  has  a 
unique  solution,  and  that  the  matrix  R = [C,  BJ  B,  0]  (first  line  : C,  B\  second 
line  : tB,  0)  of  the  system  is  symmetrical  and  sparse  with 

Cj,i  = A-f  *Pi) ) i,  I = 1,  •••,  Mh 

Bj^i  At  • <Pj  ( Si , Tj)  , j — 1,  ...,  Mh , i = 1,  *•*,  Af, 

where  for  any  u,  v € (Vh)2  , 

[u>v]om  = / £ (M)  ■ u(s,r)  ■ v(s,r)dsdr, 

Jm 

and  where  T — (a™ , ...,<*5}  , Ai,  ...,Xn)  is  the  unknown  vector.  We  obtain  a 
linear  system  RT  = L,  where  the  lines  of  L are 

/ e(M)am~1<pidsdr  + AtLam  (ipi) , 

Jm 


om  1<PMhdsdr  + (<pMh .) , 


A tai, 


Ata^. 

This  method  has  been  implemented  in  fortran,  C and  C++.  Numerical  ex- 
amples on  real  data  are  given  in  [11,18], 
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